General Disclaimer 


One or more of the Following Statements may affect this Document 


• This document has been reproduced from the best copy furnished by the 
organizational source. It is being released in the interest of making available as 
much information as possible. 


• This document may contain data, which exceeds the sheet parameters. It was 
furnished in this condition by the organizational source and is the best copy 
available. 


• This document may contain tone-on-tone or color graphs, charts and/or pictures, 
which have been reproduced in black and white. 


• This document is paginated as submitted by the original source. 


• Portions of this document are not fully legible due to the historical nature of some 
of the material. However, it is the best reproduction available from the original 
submission. 


Produced by the NASA Center for Aerospace Information (CASI) 



X-920-77-82 

PREPRINT 


A COMPUTER SOFTWARE SYSTEM 
FOR THE GENERATION OF GLOBAL 


OCEAN TIDES INCLUDING 
SELF-GRAVITATION AND 
CRUSTAL LOADING EFFECTS 

(NASA-TM-X-71308) A COMPUTRR SOfTHARB 
SYSTEM FOE THE GENERATION OF GLOBAL OCEAN 
TIDES INCLUDING SELF-GRAVITATION AND CRUSTAL 
LOADING EFFECTS (NASA) 63 p HC A04/HF A01 UucidS 

CSCL 08C G3/48 28S47 

RONALD H. ESTES 


N77-2J709 


APRIL 1977 




GODDARD SPACE FLIGHT CENTER 

GREENBELT, MARYUND 


■ 


X-920-77-82 

PREPRINT 


A COMPUTER SOFTWARE SYSTEM 
FOR THE GENERATION OF GLOBAL 
OCEAN TIDES INCLUDING SELF- 
GRAVITATION AND CRUSTAL LOADING 
EFFECTS 


R- H. Estes 

BUSINESS AND TECHNOLOGICAL SYSTEMS, INC. 


March 1977 


GODDARD SPACE FLIGHT CENTER 
GREENBELT. MARYLAND 


A COMPUTER SOFTWARE SYSTEM 
FOR THE GENERATION OF GLOBAL 
OCEAN TIDES INCLUDING SELF- 
GRAVITATION AND CRUSTAL LOADING 

EFFECTS 

R. H. Estes 

BUSINESS AND TECHNOLOGICAL SYSTEMS, INC. 

ABSTRACT 

A computer software system is described which computes global numerical 
solutions of the integro-differential Laplace Tidal Equations, Including 
dissipation terms and ocean loading and self-gravitation effects, for arbitrary 
diurnal and semidiurnal tidal constituents. The integration algorithm 
features a successive approximation scheme for the integro-differential 
system,with tliae stepping forward differences in the time variable and central 
differences in spatial variables. Solutions for M^, S^, N2, K^, K^, 
tidal constituents neglecting the effects of ocean loading and self-gravitation 
and a converged M2 solution including ocean loading and self -gravitation effects 
are presented in the form of cotidal and corange maps. 
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1.0 INTRODUCTION 


Future NASA geodynamic and oceanographic programs depend on extremely 
accurate knowledge of the geoid. Much of this knowledge will come from highly 
accurate altimeter measurements from future satellite missions ( i.e., SEASAT) 
The full use of this data necessitates realistic models of the global 
ocean tides for both the reduction of altimeter data to mean sea level and 
for direct orbital perturbations. Moreover, global ocean tide models are of 
interest in the study of energy dissipation in the earth-moon system and, for 
oceanographic studies, the global ocean tidal currents are of interest. 


The dynamics of the global oceati tides are described by the nonlinear 
Laplace Tidal Partial Differential Equations. The equations are too complex 
to solve analytically, and are amenable only to numerical solution by computer. 

In a previous work (Estes [9]) which we will refer to as report [A], a computer 
software system was developed to solve for theoretical global tides described by 
the Laplace Tidal Equations (LTE) including static earth tides, bottom friction an< 

turbulent Following Zahel [25 

technique incorporating impermeable coastal boundary conditions and a forward 
txme difference/central spatial difference numerical scheme was used. The 
present development is an extension of that work to Include ocean self-gravitation 
and crustal loading effects into the Laplace Tidal Equations. As discussed 
by Hendershott [l5l, the inclusion of these effects forces the LTE into an 
xntegro-differential system which we denote by IDLTE (Integro - Differential 
Laplace Tidal Equation) . For purposes of completene^. much of the material of 
report [A] will be repeated, particularly in the development of the tidal 

equations and in the descriptioQ of the software system design. 


Report [A] presented numerical solutions of LTE for the M, and K tidal 
constituents. We have since obtained so^l^^ S N,. K,. 0, and ^ For 

convenience «f erence, we include a complete^^^ s^^ 

solutions, which show general agreement with empirical global solutions with 
respect to the positions of amphldromic systems and tidal amplitudes . An M 
solution of IDLTE (i.e.. including ocean self-gravitat^ and crustal loading) 
has been generated and displays a marked improvement in the phase 

relations with Observations (generally wi^^^^ 30») when compared to the LTE 

2 so ution. IDLTE solutions for additional constituents will be generated 


and presented in a future publication. 

The improved comparison with observations obtained from the theoretical 
M 2 solution with the inclusion of self-gravitation and crustal loading effects 
are encouraging. However, the observed values of ocean tides contain meteorological 
effects and effects from terms beyond the second order term of the astronomical 
tidal potential which are neglected in the theoretical models. Moreover, the 
representation of the eddy viscosity deserves further study, as well as 
Internal waves, tide currents and long period ocean tides which are not 
considered in the present study. In addition, the parameters representing the 
earth model (Love numbers and load Love numbers) are averaged quantities so 
that local responses of the earth due to inhomogeneities are not included. 

These effects must be correctly modeled before agreement to the 10 centimeter 
level can be achieved. * 
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2, 0 LAPL/^CE TIDAL EQUATIONS 


The dynamic equations which describe the ocean tides, first derived by 
Laplace (1775), are 

3C 1 r ^ a *1 

3t " " Rcos^ L n cos<|.)J 

H = h(f,^) + ^ 

The first two equations are essentially a Linearized special case of the Navier- 

Stokes equations and the third is a continuity mass conservation equation. 

Here denotes latitude, A longitude, u the eastward component of velocity, 

V the northward component of velocity, and F. the components of the tidal 
force, g the gravity acceleration, R the mean earth radius, and m 
velocity of the Earth. The ocean depth h (^,, A) denotes that distance from 
the water surface to the ocean bottom when averaged over a very long time period 

■ ..'■■■ '',T; 

h«.X) . limit 

S.Cd.l.t). C^(4,l,t) denote the lnet«.t.neo„e elevatlone of the ocean .»faoe 
and oceen bottom above their time averaged , aloes. Ihos the instantaneous upward 
displacement of the ocean surface relative to the ocean bottom, the quantity 
measured by tidal gauges, is given by (see Figure 1) 


Cg((}',X,t) - C|j(<!>,X,t) 

These equations are referred to as 'long-wave" equations 
derivation the ocean is regarded as a homogeneous and incompressible fluid. 

and nonlinear acceleration terms are neglected as second order effects. iL 

assumption hydrostatic equilibrium is invote^ 

equation, in that it makes the horizontal velocity components Independent of 



the radial direction and allows integration of the expression for the balance 
of mass in a vertical water column from the ocean floor to the ocean surface. 

The vertical velocity W need not be considered explicitly , as the Coriolis term 
u)Cos(j>W is neglected in the first equation for the eastward velocity com- 
ponent. 

The displacements of the ocean and solid earth in response to the direct 
tidal forces due to the moon and sun, results in an altered mass distribution 
due to the formation of the tidal bulge inside the Earth and hence a change 
in the Earth’s gravitational field. Let and denote the gravitational 
disturbances of the bodily tide and ocean tide, respectively. Moreover,^ 
denote the forces of friction and lateral turbulent viscosity by F^ and F^. Then 

F = T + F. + F + F + F . (3) 

D S t V 

The direct tidal force, ¥, can be expressed as the gradient of the tidal 
potential T, 

T = vr W 

The gravitational disturbance due to the bodily tide results from two sources; 
a priinary response to the non— loading astronomical potential and a secondary 

response to the crustal loading by the time varying ocean tide. The Green’s 
functions for the augmented potential and elevation caused by the loading 
of the crust at an angular distance Y (spherical earth) per unit of loading 

mass are given by Farrell [ 10] as 

' ' 00 ' 

f(Y) - f S 

, e n=0 ", ■ , , . ■ 

go - . 

u’ (Y) = Aw h’ P (cosy) ' (5) 

^ ' M n n 

. ' ' ■ , e.' , ,n=0, ; ■ ■ , , 

where M is the mass of the earth and , h^ are load Love numbers. The 
primary solid earth tide, due to the fact that the free periods of the 
elastic earth vibration are less than an hour, is accurately represented as a 
static response. Thus the total body elevation and gravitational disturbance 


are 


( 6 ) 


where 


^•// 


R ?p U'(y) dn* 


\ 




C P'S*' (y) d 


(7) 


( 8 ) 


Here p is the mean density of sea-water, the domain of integration is over the 
global ocean surface and y is the angle between the integration variable 
and the position at which the integrals are being evaluated. 


COSY 


sin(j)sin(|i ’ + cos<jicos<|)* cos(X— X*). 


(9) 


The first terms (static) in the equations for expressed in 

terms of the Love numbers and hj^ with values of approximately .3 and .6, 
respectively. These parameters express the elastic yielding of the Earth 
to a non-loading potential. The gravitational disturbance of the ocean tide is 

F = 7V ^^0) 

' S , , '.S ■ ■ ■ ■ 


where 


V 

s 




GoR^gdn* 





( 11 ) 


and the domain of integration is over the ocean surface. Here G is the 
gravitational constant and x, x* represent the position on the earth s 
surface at which evaluated and the position of the element of inte- 

gration, respectively. On the surface of the spherical earth, 

Ix-x'l = JzK^H-cosy) 
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so that 


ff 


V = f I ^ d n' 

JJ 2Rsln<Y/2) 


fi 






( 12 ) 


2sin(y/2) 


In modeling the friction and lateral viscosity terms we follow Zahel [25] 
and set 


r2. 2 

\Jn 


- /«v 

f H \v/ 


F = . 

V hv 


(13) 


whete and are the coefficients of friction and horizontal edd.y viscosity > 
respectively, and A is the horizontal Laplacian 


2 2 
- 3- 1 r 

“ 2 2 2 2 2 ' 
R cos (J) 3A R 9(() 


The Laplace tidal equations thus becomes 


9 , ^ ar _g_3i 

9t Rcos(f) 9X Rcos^ 3X '^r H 


/ 2 2 

yu +v ^ o ... . 1 9r* 


u + C, Au+ 


hv Rcos(j> BA 


3t 


+ 2ojsin(!)u > 


(lHc,-h .) 




L V IL _ .s li _ 

R 3({>. R 3(j) '’r H 


C_ -4^ V + C Av + 1 


hv R 3(|) (14) 


li 

at 


. 1 r 3 (Hu) 3(Hvcos(!>) 1 _ ^ 

Rcosd. L 3X 3^ J 


V7here 


H = h(<{.,X) + ? 
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and F* is the potential due to ocean loading and self-gravitational effects. 


IP 


F’((f),A,t) = / /R C(^j)*,Ar,t)p 


Rg 1 

2sin(Y/2) 


+ (y) " 8 U' (y) 


dn' (15) 


COSY ” sin(()sin(j> ' + cos{()cos(f) ’ cos(\-A ’ ) 


The LTE thus becomes integro-diff erential in nature (IDLTE) when ocean 
loading and self-gravitation effects are included, as pointed out by Hendershott 
[L5]» Neglecting the perturbing potential F’ reduces the system to that of 
report A. The present study solves the IDLTE system by successive approximations 
using the time stepping method of Hansen Hd and finite differences, imposing 
the boundary conditions of impermeable coasts, i,e,, the perpendicular components 
of tidal velocity are required to vanish at the coasts. The procedure is to 
first generate a zeroth order solution of Equation (14) neglecting F * . Then 
F^ is calculated from Equation (15) using the zeroth order solution for the 
tide g. A first order solution of Equation (14) is then generated using the 
potential F* calculated from the zeroth order solution, etc. This successive 
iteration technique is continued until convergence is achieved. 


The periodic tidal solution is represented in the form 


* , t) = A((|) A' ) cos j^at+H^ ({()’ , 


(16) 


where contours of equal amplitude A(())',A^) are denoted corange lines and lines 
of constant phase ^ are denoted co tidal lines . The expression for F’ then becomes 


F^(4, A,t) - P((j), A) cosat - Q((j),A) sinot 


(17) 


where 


PC*)) 


.») -ffi 






sin^'d(j)'dX' (18) 
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Then 


11" _ 
dX 

ap 

n 

- ^ sinat 


11" - 
8(j) 

— cosat 

- 1 ^ 

(19) 
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3 . 0 TIDAL POTENTIAL EXPANSION 


The complex motion of the moon and sun results in a tidal disturbing force 
which displays periods ranging from thousands of years to hours. Doodson [7] 
and recently Cartwright [^] have expanded the gravitational pote^-tlal of the moon 
and sun at a point on the Earth * s surf ace into the series of tidal ” components** 
or "constituents** 


^ 6(nj^,n2,...n^)exp i(nj^T+n2S+...+ngP^) (20) 

^1*^2* ***^6 

where the sum is over all positive and negative values of the integers n^. The 
arguments of the exponential are the six astronomical frequencies presented in 
Table I. The number 


(ni ,n2+5 ,n^+5 ,n^+5,n^4-5 ,n^+5) 


( 21 ) 


is known as the Doodson numb e r an d i q uely iden t if ies individual component s . 

The functions $ are grouped into three species, for n^ = 0, 1, and 2. The 
functional forms are 

^-SPECIES;:::l^ - 


3(0, . . . ,n^) = gCKP^(sin(j)) 


SPECIES IX: 


3(1, ... = gCKP^ (sincfr) e*"^^ 


(22) 


SPECIES III: 


9 4-7i > 

6(2, = gCKP2(sln4>)e^^^^ 


The C(n„,n,,. ..n,) are constants, the p“ are Legendre polynomials and 

, 2 ■ S ; O , ■ ■ ■ , n , , . , . ■ ■ ■■ ■ 


Q M 4 

3 / moon^ 


V -J > muuiLv 

■ ^ 2 ^Me ^ „3 


- 53.7 cm 


(23) 


moon 


10 


denotes the mean lunar distance. 


where R 

moon 

The fictitious displacement 


?(t) = 


3 

g 


( 24 ) 


is denoted the **equilibrium tide*' and represents a sea surface of constant 
gravitational potential. The tidal constxtuent potential is then F — g^. 


From the definition of x, it is clear that Species III tides are semidiurnal. 
Species II diurnal, and Species I long period tides. Doodson computed the con- 
stituents with amplitude greater than lO"^ of the amplitude of the largest tide, 
SGina 400 terms. The most prominent terms are presented in Table II. 


For the present study, long period tides are not considered and for input 
purposes to the computer program, the tidal consituents are represented in the 
form . ■ 


TYPE III : 


TYPE II i! 


5 = 


a 

2 



cos (p 


gi(ot+2X) 


(25) 


5 = Ck sin2<{> 

where the constants C and a are specified by the particular constituent of 
interest,.and- ■ ■ 

0 = nj,r + n2S + n^ft + (26) 

where the rate of change is per mean solar hour. Note the factor of 2 difference 
in the semidiurnal and diurnal representation of Equation (25) . 


Definition 


Rate 

Period (radians per 

mean solar hour) 


Frequency 
(deg/solar hour) 


X Greenwich celestial 1 lunar day 
longitude 


.2529359 


Lunar orbital 
longitude 


1 tropical 
month 


.0095819 


h Solar orbital 
longitude 


1 tropical .007163 

year 


p Lunar perigee 
longitude 


8.85 Julian 
years 


.00008168 


N Lunar node 
longitude 


18.61 Julian .0000377 

years 


P Solar perigee 20,900 Julian 

® longitude years 


.344 X 10 


14.492 


.549016 


. 041068 


.004642 


-.002206 


.00000196 . 


Table I 



Designation 

Damin Doodson Number 


Amplitude Factor Frequency 
(C) (radlans/day) 


Period 

(Hours) 


Type I 


% 

(0,7,5, 5, 5,5) 

.156 

.459969 

327.84 

Type II 

0, 

X . ■ 

(1,4,5, 5, 5,5) 

.377 

5.840445 

25.82 


(1,6, 3, 5, 5, 5) 

.176 

6.265983 

24.07 

* 4 . 

(1,6,5,5,5,5) 

.531 

6.300388 

23.93 

Type III 

«2 

(2,4,5. 6, 5, 5) 

.174 

11.910864 

12.66 

M2 

(2, 5, 5, 5, 5, 5) 

.908 

12.140833 

12.42 

^2 

(2, 7, 3,5, 5, 5) 

.423 

12.56637061 

12.00 

•^2 

(2, 7,5,5, 5, 5) 

.115 

12.600776 

11.97 


Table II 


A . 0 COI-fPlTTATIONAL TECHNIQUES 


4 . 1 Haasen Grid 

The numerical technique employed in the present study for the Laplace Tidal 
Equations (14) follows Hansen [12] and Zahel [25] in applying a finite difference 
algorithm in temporal and spatial variables and integrating forward in time 
( timestepping) . We employ the global grid system proposed by Hansen (also known 
as the Richardson grid) whereby a staggered grid with velocity components and 
tidal elevations are computed at adjacent, but not coincident, grid points. The 
coastal boundary conditions are represented by horizontal and vertical lines 
imposed on the grid and pass only through those points denoting velocity com- 
ponents. Vertical boundary lines pass through u grid points (eastward velocity) 
and horizontal boundary lines pass through v grid points (northward velocity) so 
that the requirement of the vanishing .of perpendicular components of velocity 
forces, all velocity components falling on the boundary line to be Identically 
zero for all time. Moreover, it is to be noted that with this procedure there 

is never a corner point involved in the boundary conditions. This is of consid- 
erable numerical importance. The Hansen grid is depicted in Figure 2. 


4.2 Finite Difference Equations 

The finite difference formulation for the Laplace Tidal Equations (14) , using 
central differences for spatial variables and forward differences for the time 
variable are. following Zahel [25] 


. 2 , 1 


1 C _At ^/l(N ,M, t_) -^[ v(N ,t„)+v(N ,Mvt)+v(N+l ,H-1 )+v(N+l,M,t ) ]/ 

~ ^ n . n n 

I |[h(N,M)+h(N,M-l)+C(N,M,t^)+5(N,M-l,t^)] 


• C^) + 2uAt sin[!{i(N) ] 


i [ v(N , M-1 , y +v(N , M, t^) +y (N+1 , 


+ C. At 
hv 


( u(N,K-l,t )+u(N,M+l,t )-2u(N,M,t ) 
^ n n n 

V [Ilcos[4i(N)lAX]^ 


u (N+i ,M, t )+u(xN-l, M,t )-2u(N ,M, t ) 
n n 


[RA<},]' 


■"! 


14 


Rcos [ (N) ] 




jf(N,M)cos[at^] - f (N,M)sin[atJ j 


+ F^ht 


j C^At ^/^N,M.t^)^+ ^[u(N-l , M, t^)+u(N-l .^H-l, t^)+u(N,m, t^)+u(N,M,t^) ! 

i[h(N,M)+h(N-l,M)+ 5 (N,M/^^^ 


, v(N,M, t^) - 2 ojAl: sin[ 4 >(N) - -“1 
. A[u(N-X,M,t )+u(N-l,m,t_)+u(N\>Hlvt^^^ 

tf n 


( v(N,M-l,r )+v(N.XH-l,tJ -2v(N,M.t^) (27 

( [Rcos it(N)^ 

:n +1 ,M, t^)+v(K-l,II, t^)- 2 v<N ,M, t^) V 

' / : : 


ML 

^ I 




SP _ P(N.M) - f (N.M-1). 

— (N.M) - y 


(N.M) * 


)(N.M) - Q(N.M-l) 
AX 


(N,M) = 


P(N.MV - P(N-l.M) 


(N,M) - 


15 


where 


■p s: S 

A Rcos^ 8A 


p - £ IL. 
<p R S<t) 


and 


= 


CSK 


Als 


^ cos tMN)] sin [2(A(M) - ^) + ot^] 


SPECIES IIIJ 


^<P 


sin [2((1)(N) - -^)] cos [2A(M) + ot^] 


sin [(|,(N)] sin [(A(M) - % + ot ] 


(29) 


SPECIES II: 


F. = cos [2(<j.(N) - -^)] cos [A(M) + at^] 


gp 9Q 3P 90 

The quantities , and arise from the perturbing potential F\ and 

are held constant in time over an iteration (see Section 2.0). For computing efficiency 
these derivatives are precalculated and held in storage for the integration 
algorithm. The continuity equation becomes 


At , f th(N,M)d*h(N,M~l)+e(N,M;t^)+C 
^ Rcos [(j)(N) ] I 


AA 


f[h(N,m'l)‘hh(N,MH^ 


AA 


~[h(N;M)+h(N-l,^0^-e(N,M;t )+C(N-l,M,t^)]v^^^ 


AA/cost^(N)“- 


|[h(N+l»M)+h(N,M)+C(N+l,M, t^)+5(K*M, t^) 3v(N+l,M, j 


A^/cos [ «j) (N+1) - ] 

16 


( 30 ) 



/ 



hH M X(M) M4-1 




Figure 2 
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is calculated in 


It is to be noted that the tidal elevation C at time t^^^^ 
terms of the velocities at time finite difference formulas for the 

velocities u and v involve only nearest neighbor points ^ with the exception of 
the Laplacian terms for turbulent viscosity. For grid points adjacent to coastal 
boundaries, the evaluation of the Laplacian central difference term involves a 
grid point which would fall on land. The computer algorithm handles this by 
assuming the normal derivative of velocity at the boundary to vanish, thereby 
eliminating the need to evaluate the point. In this way values at land grid 
points are never involved in the algorithm. (Coastal boundary points are not 

considered to be land points.) Note that these Laplacian terms are small relative 

2 . ■ ■ ■ 

to the other terms, having an R divisor . 

4.3 Stability 

The stability of the system desqribed by Equation (lA) neglecting 
the r’ terms has been discussed in Report A. In particular, it was shown 
that the global tides computed by the time stepping procedure become periodic 
in time and independent of the initial state (providing the initial values 
of u, v, 5 were small enough that the system did not iimnediately blow up) 
within five tidal periods, provided the time step At is small enough to 

insure stability. The strong dissipation of the friction and lateral 

viscosity terms are the cause of this rapid convergence, quickly damping 
out the initial transients in the numerical integration and eliminating 
spurious waves. A stability analysis for Equations a4> less the rV terms 
has been performed by Zahel EZ5l where formulas for the maximum time step 
At as a function of grid spacing, latitude limits and dissipation parameters 

C and C, are given which will guarantee stability for the tidal solution. 

r hv . . ^ 

Report A describes a numerical procedure for easily establishing a suxtable 

time step. ■ 

The successive approximation method proposed in Section 2.0 presents 
no difficulty in establishing a for a given iteration to insure stability, 

as the r’ terms on the right hand side of Equations are known functions of 
position and time. The conditions under which a successive approximation 
procedure will converge to the solution of the integro-differential system 
described by Equations 0.4) for the perturbing potential T • will not be 
pursued in this report. 
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4.4 totegratlon Algorithm and Solution Representation 


The program implementation of the integration algorithm for Equations 

(14) and their boundary conditions is briefly described in Report , A. The basic 

design is such that calculations are performed on a single latitude band at 

a time and required quantities are accumulated over three bands. The computer 

code features compact arrays for all variables u, v, g and h with values 

only at non-boundary ocean points on the global grid. Compact arrays also exist 
9P 9P 80 30 

f or — , ^ precalculation of the derivatives of P and Q Increases 

the core requirements of the program, but results in a substantial savings of 
CPU time. All variables required by the integration algorithm reside in core 
which helps keep I/O and CPU time requirements down. The time needed to integrate 
Equations on a global 3® x 3® grid with a 3 minute step size for 12 hours 
is approximately 3 CPU minutes on an IBM 360/91. Thus, to integrate over 6 
periods of a semi-diurnal tide requires approximately 20 CPU minutes. The CPU 
time increases substantially for finer global grid meshes (and consequently 
smaller time steps) • 


The program is designed with an overlay structure to minimize computer core 
storage requirements. A solution with a global grid resolution of 3° x 3® 
can be accommodated in under 500K bytes of core. A 2® x 2® global grid requires 
approximately 600K bytes of core, while a global 1® x 1® grid would need 
slightly in excess of lOOOK bytes. 

The tidal solution ?(<|),X;t) obtained by the software system in time-stepping 
from time t=0 to t=T is specified on the global grid at time t = T. It 
is convenient (and conventional) to represent the periodic solution in the 
form (Equation (16)) 


?(4>,X;t) = A((|),X) cos [>!'((f,X) + at] 
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where t is expressed in mean solar time. The computer software establishes 
this representation by integrating the solution over a quarter of a tidal period 
from t = T to t - T + V2o so that 


A((j.,X) = J e(f.X5T)^ + ?(<!>, A;T+TT/2a)2 


= tan" 


-£((fc.X;T+ /2q 


- 0T 


(31) 


The software system employs a graphics package to display solutions in 
terms of corange lines (contours of equal amplitude A(<{),X)) and co tidal lines 
(lines of equal phase 'F((j>,X)). For purposes of display, we employ the representation 


— A((j) , X) cos [ot“S ((() , X) ] (32) 

where 

S((j),X) - - 'F((|),X) , (33) 

The cotidal values appearing on the tidal charts presented in this report are 
expressed in hours 0 to 24, instead of angular measurement. Hour values are 
obtained by dividing the phase function S expressed in degrees by 15 degrees per 
solar hour. Note that both diurnal and semi-diurnal tides have cotidal values 
expressed in hours from 0 to 24. These relationships are given explicitly as 
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4.5 Rvaluation of Green's Functions fo r Ocean Loading. 

The Green's functions for the augmented potential and elevation caused 
by the loading of the crust given by Equations (5) are slowly converging 
series, as discussed by Farrell M . Following Farrell, and using the 
identities 

g P^(COSY) = -^(y72) 

= - £n|2sin(Y/2) jxl+sin(Y/2)j > 

^ n \ ' 


we rewrite Equations (5) in the form 


"■(Y) - ir { £ 'n * '•« S'" 


X I + V (h'-h„ ) p_ (cosy) i 

Mel2sin(Y/2) S “ ) 


(35) 


$» 


(Y) - f { S C*° n~ -) 

e \ n-0 


+ k 


P (cosy) ) 

"^1 




j ^ "^n ~ ^ P^(cosy) - k Zn r2sin(Y/2) [1 + sin(Y/2)j 


where 


b = lira h’ 
"oo n-^ n 


Xp * 4^ ^n 


( 36 ) 


We have investigated several summation techniques for evaluating the series 
in Equation (35), including the Shanks coirvergence acceleration method £23]. 
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However, we found the Euler transformation method discussed by Farrell to be 
the most efficient. The values of and h^ used are those given by Farrell 
for the G-B Earth Model, where 

k = - 2.482 


h = - 5.005 


Moreover, we have used the value 


ho = - .134 


obtained by Longman M tor a Gutenberg Earth Model. Farrell gives values of 


h’, k^ for 
n n 


n = 1 , 2 , 3 , 4 , 5 , 6 , 8 , 10 , 18 , 32 , 56 , 100 , 180 , 325 , 550 , 1000 , 1800 , 3000 , 10000 . 

Quadratic polynominals were fitted through these values to obtain k^, h; for 
other values of n required for the sximmations . 

Table III presents calculated values for U'(y). $’(y) and the Greens function 


R^P 




(37) 


az 


For convenience of comparison, the values of U* and are scaled by (Ry) x 10 , 
with Y in radians. The quantities are in mks units, and an applied load of 
1 kg. is assumed. Values taken for R, P, g and M^ are 6371000., 1030., 9.82, 
and 5.975 x 10^\ respectively. 
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TABLE III 


Green’s Eunction 



1” 

-13.196 

-1.663 

73.178 

2" • 

- 9.725 

-1.159 

30.277 

3° 

- 7.490 

-0.025 

17.576 

4“ 

- 6.028 

1.368 

11.964 

5“ 

- 5.294 

2.886 

9.144 

6“ 

- 4.625 

4.494 

7.310 

r 

- 4.164 

6.158 

6.112 

8® 

- 3.971 

7.850 

5.339 

•:9® : : 

- 3.921 

9.580 

4.798 

O 

o 

- 3.718 

11.367 

4.311 

0 

O 

- 2.603 

31.214 

2.328 

30® 

- 1.577 

53-987 

1.716 

o 

O 

< 

- 0.293 

78.671 

1.406 

Ln 

O 

0 

0.838 

104.017 

1.238 

60® 

ll674 

128.809 

1-142 

O 

6 

2.076 

152.136 

1.089 

O 

O 

00 

2.047 

173.445 

1.061 

VO 

o 

o 

1.647 

192.523 

1.046 

100® 

0.924 

209.427 

1.039 

no® 

- 0.029 

224.447 

1.035 

120® 

- 1.138 

238.011 

1.033 

130® 

- 2.311 

250.661 

1.032 

140® 

- 3.439 

262.991 

1.030 

150® 

- 4.506 

275.579 

1.028 

160® 

- 5.438 

288.956 

1.027 

170® 

- 6.263 

303.668 

1.027 

180® 

- 6.660 

320.324 

1.025 



4.6 Tidal Solutions Neglecting Loading and Self-Gravitation Effects 


Solutions for the principal tidal constituents K 2 , K^, 0^ 

and have been generated using the integration algorithm described in 
Section 4.1 neglecting the terms involving F’ and are displayed in Figures 3-9. 
These tide models are then solutions of the LTE as described in Report A, 
and do not exhibit an integro-dif ferential character (IDLTE) . They would, 
however, serve as the Initial approximations for the successive approximation 
algorithm described in Section 2.0 for the IDTLE (see Section 4.7). Global 
grid resolutions of both 3° x 3® and 2® x 2® have been used for generating 
tide models for the principal constituents and display excellent agreement. The 
Love parameter values used in these solutions were = .3 and h^^ - .6 while 
the dissipation parameters were “ .003 and ~ 10^, respectively (the 
same values used by Zahel [25, 26]) • In the 2® x 2® numerical models the 
effects of the tides of polar seas north of 80® latitude are ignored 
(81® latitude for the 3® x 3® models) . This is due to the fact that increasing 
the northern latitude licoit necessitates decreasing the integration step 
size to maintain stability, hence greatly increasing the computer time needed 
to arrive at the solution. However, our numerical experimentation has shown 
that increasing the northern latitude limit has no discernable impact on the 
solution south of 80®. The ocean depth profile used in the solutuions is 
presented in Figure lOi Initial values of u, v and^ at time t = 0 were set 
to zero. A time series plot of the tidal solutuion, g, for a grid point 

in the central Atlantic is presented in Figure 11 . The solution grows from 
the zero initial solution to a steady state solution within several periods. 

All solutions for the principal constituents presented in this section were 
integrated over eight tidal periods (assuring a steady state solution) before the 
amplitudes and phases (Equation (31) ) were calculated. 


Prominent features evident in both the diurnal and semi-diurnal solutions 
are amphidromic systems, or geographic regions where the tidal amplitude 
vanishes. Cotidal lines radiate from the amphidromic points, where the sense of 
the amphidromic system is denoted as the direction in which the high tide 
occurs progressively later . Amphidromes are usually (not always) counter- 
clockwise in the Northern Hemisphere and clockwise in the Southern Hemisphere. 
These stationary centers can result physically as a consequence of the character- 
istic wave velocity being of the same order of magnitude of the velocity 
the sublunar point. 
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The family of semi-diurnal waves S^, N^, exhibits very similar character- 
istics in the tidal solutionSj> as does the family ^ 1 * ^1 diurnal waves* 

For this reason we will only discuss the and solutions of Figures 3 and 
7 in detail. 

The M 2 solution of Figure 3 was compared with the numerical solutions of 
Hendershott [151, Pekeris and Accad EH, Zahel E5], and Bogdanov and Magarik [ 3 ] as 
to tidal amplitude, phases and location of amphidromic systems .These solutions are con- 
veniently summarized by Hendershott [14], Our M 2 solution compares very closely 
With Zahel, with nearly a systematic 30% decrease in tidal amplitudes consistent 
V 4 -th the factor (l+k 2 -h 2 ) in our solution. The various models are in fairly good 
agreement in the North and South Atlantic, with the exception that Hendershott 
does not find the South Atlantic araphidrome, and the atnphldrome at the southern 
tip of South America is farther west in ours and Zahel* s solution relative to the 
dissipative model of Pekeris. All models are again quite similar in the Indian 
Ocean, being dominated by a large **antiamphidrome*^ i, e. a region of high tidal 
amplitude moving nearly synchronously. The amplitude of this feature shown by 
Hendershott is substantially larger than the other models. Our experimentation- 
SiiOws that the region at the southern tip of Madagascar is near an amphidromic 
state. The Pacific basin shows the greatest differences between the numerical 
models. Our solution, along with Zahel, shows the appearance of an amphidromic 
s>suem off of California, and a line of rapid phase change in a concave arc 
running from Central America to Southern Peru. Our model, in agreement with Zahel, 
shcistwo antiamphidromes on either side of the dominant South Central Pacl fie 
araphidrome, although the eastern antlamphidrome of our model is more pronounced 
than that of Zahel. The Western Pacific .basin does -not show the pseudo-nodal 
line extending from Japan to New Guinea present In the empirical chart o£ 

Dietrich [6]. This region in our M 2 model Is instead dominated by an amphi- 
dromic system. The full 360“ counter clockwise phase change about New Zealand 
is in agreement with the empirical charts. 

The Ki tidal solution in Figure 7 was compared with that of Zahel [2^. 

The amphidromic Systems of the Atlantic and Indian Oceans correspond quite closely 
between the two solutions, with the exception that the zero cotldal contour 
runs directly between the North Atlantic and Central Atlantic amphidromes of 
our model rather than to the coasts of Spain and Dakar, respectively. As with the 

^2 hhe Pacific basin offers the most contrast between the numerical models . 
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Where Zahel finds a single dominate amphidrome in the Central Pacific, we find 
this region divided into two centers along an arc running from Wake Island 
to Tahiti. In addition we find an amphidromic system off the coast of Chile, 
two models compare closely in the northeast Pacific and the South Pacific. 

As with the wave, general agreement is found with the empirical chart of 
Dietrich with the exception of the western and southeastern Pacific basin. 

In particular, the region of constant phase in the western Pacific found by 
Dietrich is dominated by an amphidromic system in pur model. 

Tables IV-V compare tidal amplitudes and phases of the 0^ 

models at selected locations with the harmonic constants available from the 
Admiralty Tide Tables [1]. The phases are relative to the Greenwich meridian, 
and are the S((j), y) phases defined by Equation (3$ . 



The 
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Figure 3 

2* X 2* K 2 (Neglecting Loading and Self-Gravltatlon Effects) 


S; Dashed lines cotldal lines (phase) in hours 

A: Solid lines corange lines (amplitude) in decimeters 

5(4., X;t) - A(4>,X) cos lot - S] 

0 - .00014052 radians/sec 
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Figure 4 

2® X 2® Sp (Neglecting Loading and Self-Gravitation Effects) 

S: Dashed lines cotidal lines (phase) in hours 

A: Solid lines corange lines (amplitude) in decimeters 

C(4>,^;t) * A((|),X) cos[ at - S] 

o * ,0001454452 radians/sec 
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Figure 5 

2® X 2® N 2 (Neglecting Loading and Self-Gravltatlon Effects) 

S: Dashed lines co tidal lines (phase) In hours 

A: Solid lines corange lines (amplitude) In centimeters 

£(4>»A;t) ■ A(4>,A) cos [at - S] 

a • .00013788 radlans/sec 





Figure 6 

2* X 2* K 2 (Neglecting Loading and Self-Gravitation Effects) 

S: Dashed lines cotidal lines (phase) in hours 

A: Solid lines corange lines (amplitude) in centimeters 

C(().,X;t) “ A(4(,X) cos [ot - S] 

o “ .00014584 radians/sec 
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Figure 7 

2* X 2® (Neglecting Self-Gravitation and Ocean Loading Effects) 

S: Dashed lines cotidal lines (phase) in hours 

A: Solid lines corange lines (amplitude) in decimeters 

C(<|>,A;t) * A(<^,A) cos [ot - S] 

o » .0000729216 radians/sec 



Figure 8 

2® X 2® 0^ (Neglecting Self-Gravitation and Ocean Loading Effects) 

S: Dashed lines cotidal lines (phase) in hours 

A: Solid lines corange lines (amplitude) in decimeters 

t(4>,X;t) = A((t>,X) cos [at - S) 

o * ,0000675983 radians/sec 






Figure 9 

2® X 2® (Neglecting Self-Gravitation and Ocean Loading Effects) 

S: Dashed lines cotidal lines (phase) in hours 

A: Solid lines corange lines (amplitude) in decimeters 

C(<t>.A;t) » A(<|>,X) cos [at - S] 

a “ .0000725236 radlans/sec 












TABLE IV 


(Ti 


Comparison of and Tide Models with Observations 
(Neglecting Loading and Self-Gravitation Effects) 



OBSERVED 

M 2 MODEL 

S^ OBSERVED 

S^ MODEL 

Azores (N. Atlantic) 

63.3^ 

26.3* 

82®. 5 

48® 

(37* 44* N. 25* 40’ W) 

49 cm. 

39 

18 cm. 

18 

Bermuda (N. Atlantic) 

20^ 

300* 

45® 

310® 

(32* 19* N. 64* 50' W) 

38cm. 

19 

8 cm. 

8 

St. Helena (S. Atlantic) 

87® 

54* 


70® 

(15* 55’ S. 5* 42* W) 

34 cm. 

25.5 


11 

Chagos Is. (Indian Ocean) 

270® 

252“ 

300® 

268® 

(7* 21’ S. 72* 27’ E) 

52 cm. 

26 

27 cm. 

13 

Cocos Is. (Indian Ocean) 

306® 

266° 

5® 

288® 

(12* 5’ S. 96* 51’ E) 

27 cm. 

22 

9 cm. 

12 

San Francisco (E. Pacific) 

215® 

195" 

210® 

200® 

(37* 48’ N. 122* 27’ W) 

54 cm. 

35 cm. 

12 cm. 

14 

Honolulu (Pacific) 

61® 

35 

62® 

63® 

(21* 18’ N. 157* 52’ W) 

16 cm. 

33 

5 cm. 

13 

Easter Is. ( S. Pacific) 

32® 

281* 

240® 

300® 

(27* 9’ S. 109* 27’ W) 

20 cm. 

35 

6 cm. 

16 

Apia, Samoa (S. Pacific) 

171® 

87* 

150® 

105® 

(13* 48’ S. 171* 46’ W) 

38 cm. 

54 

9 cm. 

24 

Siapan (S . Pacific) 

300® 

159* 

315® 

190® 

(15* 12’ N. 145* 43’ E) 

17 cm. 

21 

5 cm. 

9 

Yaruto, Marshall Is. (S. 
Pacific) 

125® 

88* 

150® 

109® 

(5* 55’ N. 169* 39’ E) 

47 cm. 

46 

26 cm. 

23 

Burgana, Solomon Is. (S. 
Pacific) 

134® 

196* 

180® 

210® 

(9* 11’ S. 160* 13’ E) 

11 cm. 

33 

8 cm. 

15 



TABLE V 


Comparison of and 0^ Tide Models With Observations 
(Neglecting Loading and Self-Gravitation Effects) 



OBSERVED 

MODEL 

0^ OBSERVED 

0^ MODEL 

Azores (N. Atlantic) 

66 ® 

70® 

315® 

80® 

(37® 44^ N- 25*^ 40* W) 

4 cm. 

4.6 

3 cm. 

2 

Bermuda (N. Atlantic) 

195 

230® 

210 ® 

220 ® 

(32® 19* N. 64® 50* W) 

7 cm. 

4 

5 cm. 

3.5 

S t . Helena (S , Atlantic) 
(15® 55* S. 5® 42* W) 

■ , 

165® 

3 

- 

137® 

2 

Chagos Is, (Indian Ocean) 

84® 

80® 

60® 

62® 

(7® 21’ S. 72® 27* E) 

3 cm. 

4 

3 cm. 

3 

Cocos Is • (Indian Ocean) 

165® 

168® 

127® 

139® 

(12® 5* S, 96® 51* E) 

11 cm. 

11 

9 cm. 

8.5 

San Francisco (E, Pacific) 

228® 

229® 

210 ® 

193® 

(37® 48* N. 122® 27* W) 

37 cm. 

43 

23 cm. 

29 

Honolulu (Pacific) 

225® 

266® 

225® 

230® 

(21® 18* N. 157® 52* W) 

15 cm. 

19 

8 cm. 

15 

Easter Is, ( S . Pacific) 

18® 

57® 

90 

17.5® 

(27® 9* S. 109® 27’ W) 

8 cm. 

8 

5 cm. 

8 

Apia, Samoa (S . Pacific) 

75® 

28® 

60® 

340® 

(13® 48* S. 171® 46’ W) 

5 cm. 

7 

3 cm. 

3 

Siapan (S. Pacific) 

75® 

130® 

45® 

82® 

(15® 12’ N. 145® 43’ E) 

14 cm. 

11 

10 cm. 

11 

Yaruto, Marshall Is. (S. 

75® 

306® 

45 ® 


Pacific) 

(5° 55* N. 169° 39' E) 

9 cm. 

3 

6 cm. 

2 

Burgana, Solomon Is. (S. 

45° 

60® 

22 ® 

57® 

Pacific) 

(9° 11’ S. 160° 13’ E) 

20 cm. 

8 

11 cm. 

6 


4.7 


Tide Including Ocean Loading and Self-Gravitation 


The SGlutions of Section 4.6 show general agreement with empirical solutions 
with respect to the positions of amphidromic systems and tidal amplitudes. 

Notable exceptions to the agreement of the theoretical with empirical solutions 
are in the western and southeastern portions of the Pacific basin. Phase relations 
also differ in portions of the North and South Atlantic. To improve the 
model of Section 4.6, the method of successive approximations outlined in 
Section 2.0 was used to solve the IDLTE given by Equation 0-4) on a 3® x 3° 
global grid. The solution of Section 4.6 was used as the initial solution. 

The convergence of the procedure was slow in the southeastern Pacific 
basin, requiring 16 iterations to reach final convergence. Convergence in all 
other regions was more rapid, requiring approximately 8 iterations, (con- * 
vergenc^^: was considered to be achieved when the corange and co tidal lines 
no longer changed from one iteration to another) . To assure that the iteration 
procedure had stabilized, two additional iterations (18 iterations) were 
calcualted. Figure 12 displays the self- gravitation and ocean loading 
potential P’ obtained by evaluating Equations (L5) with the initial solution 
from Section 4.6 ( the Gth order approximation). The portions of P' due to 
self-gravitation only and crustal deformation only (U* Green’s Function) are 
displayed in Figures 13 and 14 respectively . The converged solution (16th 
iteration) is presented in Figure 15 . 

Comparison shows no large differences in amplitude or position of amphidromes 
between the solution of Section 4.6 and that of Figure 15 which includes 
the perturbing effects of ocean self ^gravitation and crustal loading. There are, 
however, significant changes in phase. In particular, the phases of the North 
Atlantic amphidrome are rotated clockwise approximately 25® and the phases 
from the South Atlantic amphidrome are rotated clockwise by 45® along the 
South American coast and approximately 20® along the equatorial coast of 
Africa. These shifts in phase improve agreement with observations in the 
Atlantic as demonstrated by the comparisons at Bermuda, the Azores and St. 

Helena in Table VI • The rotation of phas es by approximately 20 ® in a counter 
do ckwi s e manner about the amphi dr ome off the wes t c o as t o f Aus t r alia 1 ikewi s e 
improves agreement with observations in the Indian Ocean , as demonstrated by 
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Figure 14 

3® X 3® Crustal Response to Ocean Tidal Loading (ON ZEROTH ITERATION) 

S: Dashed lines phase lines in hours 

A: Solid lines amplitude lines in centimeters 
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Figure 15 

3* X 3® (Including Loading and Self-Gravitation Effects) 

S: Dashed lines cotidal lines (phase) in hours 

A: Solid lines corange lines (amplitude) in decimeters 

C(*,A;t) =* A(({),X) cos [ot - S] 

o * .00014052 radians/sec 
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K, OCEAN TIDE PAHSE AND AMPLITUDE 



OBSERVED 

THEORETICAL 

THEORETICAL 



(Without Ocean Loading 

(With Ocean Loading 



and Self-Gravitation) 

and Self-Gravitation) 

Azores (N. Atlantic) 

63. 3*^ 

26.3“ 

50° 

(37° 44' N. 25° 40’ W) 

49 cm. 

39 

41 

Bermuda (N. Atlantic) 

20^ 

300° 

308° 

(32° 19’ N. 64° 50' W) 

38 cm. 

19 

17 

St. Helena (S. Atlantic) 

87® 

54° 

73° 

(15° 55’ S. 5° 42’ W) 

34 cm. 

25.5 

26.2 

Chagos Is. (Indian Ocean) 

270® 

252° 

269° 

(7° 21’ S. 72° 27’ E) 

52 cm. 

26 

31 

Cocos Is. (Indian Ocean) 

306® 

266° 

285° 

(12° 5’ S. 96° 51’ E) 

27 cm. 

22 

29 

San Francisco (E. Pacific) 

215® 

195° 

210° 

(37° 48’ N. 122° 27’ W) 

54 cm. 

35 cm. 

45 cm. 

Honolulu (Pacific) 

61® 

35 

66° 

(21° 18’ k. 157° 52’ W) 

16 cm. 

33 

31 

Easter Is. (S. Pacific) 

32® 

281° 

303° 

(27° 9’ S. 109° 27’ W) 

20 cm. 

35 

34 

Apia, Samoa (S. Pacific) 

171® 

87° 

112° 

(13° 48' S. 171° 46’ W) 

38 cm. 

54 

55 

Slapan (S. Pacific) 

300® 

159° 

200° 

(15° 12’ N. 145° 43’ E) 

17 cm . 

21 

25 

Yaruto, Marshall Is. (S. Pacific) 

125® 

88° 

114° 

(5° 55’ N. 169° 39’ E) 

47 cm. 

46 

53 

Burgana, Solomon Is. (S. Pacific) 

134® 

196° 

210° 

(9° 11’ S. 160° 13’ E) 

11 cm. 

33 

34 







comparisons at Chagos Island and Cocos Island. In the Pacific northeast, 
the phase lines radiating from the amphidrorae off the California coast are 
rotated clockwise by approximately 15® improving the model along the North 
American western seaboard. The only appreciable shift of an amphidrome involves 
the system off the coast of Equador, which moves southward by approximately 
15® of latitude as well as being rotated clockwise by 45® . The amphidrome which 
dominates the northwest Pacific basin experiences a 45® counter clockwise 
shift of the phase lines along the Asian coast making phases along Japan 
and Okinawa more realistic. Agreement of the M 2 solution including ocean 
loading and self- gravitation effects with island observations is generally 
within 30® of phase and 25 cm. is amplitude. Exceptions are along the Aleutions, 
in the Sargasso Sea region and the southeastern Pacific basin, (see Table VI). 

Solutions for and S 2 including the effects of ocean loading and self- 
gravitation are currently being calculated and will be presented in a later 
publication. 

4.8 Total Tide Comparison with Observations 

The complexity of the total tide becomes evident when the seven constituent 
models are combined to yield an approximation to the total ocean tide. The 
theoretical global ocean tide models are not expected to be realistic along 
highly irregular coast lines or in shallow water along continental coasts due 
to the assumptions imposed in the boundary conditions. Deep ocean islands, 
such as the principal ports of Honolulu, Bermuda and the Azores of the Admiralty 
Tide Tables provide a reasonable means of comparing the ocean tide models 
with observations. It must be remarked, however, that groups of islands, 
to small to be included in the global model, can strongly influence the tides in 
their vicinity. Figures 16-19 display the tide composed of the S 2 , 

^2* ^1’ ^1* ^1 Section 4.6 and the M 2 model of Section 4.7 for 

the Azores, Bermuda, San Francisco and Honolulu. Times and values of high 

and low water from the Admiralty Tide Tables are indicated on the plots. 

The comparison at the Azores, which is dominated by the M^ Tide, is very 
good, while the predicted tide at Bermuda displays a reasonable phase but 
small amplitude (see also Table VI) . The tides at San Francisco and Honolulu 
show a greater diurnal influence. The phases and amplitudes of the predicted 
tide at San Francisco are reasonable, while the predicted tide at Honolulu shows 
poorer agreement. Causes for possible disagreement between observations and 
the theoretical models were discussed in Section 1.0. 
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Figure 16 


Tide Model Prediction for Azores, Jan. 1, 2 1977 
(Circled Points are Observed High and Low Waters) 



Hours, GMT 


Meters 


Figure 17 

Tide Model Prediction for Bermuda, Jan. 2 1977 
(Circled Points are Observed High and Low Waters) 



Hours GMT 


Figure 18 

Tide Model Prediction for San Francisco, Jan. 1, 2 1977 
(Circled Points are Observed High and Low Waters) 



Hours, GMT 






Figure 19 

Tide Model Prediction for Honolulu, Jan. 1, 2 1969 

(Circled Points are Observed High and Low Waters) 



Hours j GMT 


5.0 Computer Software System 


The integro-differential IDLTE Software System is composed of three 
basic modules; 

I. IDLTE Integration Package. This program itegrates the system 
described by Equation (14) of Section 2.0. It requires as 
input the prescribed continental boundaries, the global 

ocean depth profile defined on a grid, an initial solution, and 
the potential TV. To calcualte the zeroth order solution both 
the initial solution and F’ are set to zero. 

II. Surface Integration Package. This program evaluates the potential 
r* of Equation (15) of Section 2.0. It requires as input the 
prescribed continental boundaries, a specified global tide 

5((t),X) in the form of gridded amplitude and phase values and a 
table of values for the Green’s functions^’ and U’ of Equation 
(5) as a function of angle . 

III. Ocean Loading Green’s Function Package. This program evaluates 

a table of the Green’s functions and U’ of Equation (35) from 

an input set of load Love numbers k’ and h’ . 

n n 

The three modules are written in FORTRAN IV, and Modules I and II 
are designed with a modified graphics package (WOLF PLOT PACKAGE) for solution 
display. The I/O design of the system uses multiple magnetic tape drives. 
Table VIII defines the I/O units for Modules I and II. 

5.1 Program Input 

User control of Modules I and II is set up by card input for the 
princij- 1 yariables through a NAMELIST specification. The variable names, 
their function and their default values are presented in Tables IXand X . The 
ocean depth values required for Module I are read from an external tape, with 
values defined on a world grid (depth values on ocean grid points and 
0 on land grid points) . 

Following the NAMELIST INPUT for both Modules I and II, the continental 
boundary data is input . The ocean boundaries are ins ti tuted as horizontal and 
vertical segments on the grid with horizontal boundaries along v velocity 
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points and vertical boundaries along u velocity points. The coastal boundary 
configuration is input by cards in a special compact format. For each 
latitude on the world integration grid there is a data card of the format 
1415 with fields. 


LAT NO IDATA(12) 

where LAT is the latitude in degrees, NO is an integer 0 or 1 depending on 
whether the point at zero longitude at the particular latitude is on land or 
ocean, respectively, and the array IDATA is the value of the longitude, in 
degrees, at each grid point for that latitude where a land~weter boundary 
occurs. With this input technique any arbitrary coastal boundary configuration 
may be implemented. The boundary configuration presently employed is a 
realistic 1° x 1° world ocean boundary. The user specified grid size extracts 
the proper ocean boundary from this I" x i« base. 

For Module II, the ocean boundary data is followed by a table of values 

for the Green's Functions U' (y) and$' (Y) for y=o'*,1®, 2®, . .., 180® 

—12 ' 

(scaled by (Ry) x 10 as discussed in Section 4<5). The input format is 
8F10.3. 

Module III is a stand-alone program designed to calculate the Green's 

Functions and U' as described in Section 4. 5 . The program need only be 

run when a new set of load Love number h' and k' are employed to calculate 

n n 

new functions and U' . The data for the program consists entirely of values 
and h^, for a set of values of n. This data is input in data 
statements within the main program of Module III, 



Subroutine Name 
SPLASH 


MAINA 


TI'DPLT 

FUNS 

ALGR 

DEPTH 

STOITE 

ZARY 


Table VII 

MODULE I SUBROUTINES 

Function 


The executive program which reads Namelist Input and calls the 
working subroutines of the system. This subroutine contains 
a BLOCK DATA subroutine. 

The module which reads the input coastal boundary data and sets 
up the world function (0 on land points and 1 on ocean points) 
corresponding to the selected integration grid size. Prints 
the world function to display boundaries and defines the trans- 
formations from the global grid to the compacted ocean arrays. 

The plotting module which generates corange contours, cotidal 
lines and time series plots for selected geographic points. 

The function subroutine which evaluates the derivatives of the 
tidal component for use in the driving terms of the Laplace 
Tidal Equations. 

The integration algorithm for the Laplace Tidal Equations. 

The module which transforms the input ocean depths defined on 
a world grid to the compacted ocean grid. 

INPUT /OUTPUT routine for the world grid ocean depth values. 

The module which transfonas the compacted ocean arrays onto 
world grids for display and plotting purposes. 
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Table VIII 

TAPE I/O UNITS FOR MODULES I & II 


Module I/O Unit 


Function 


Module I 30 

40 
50 

60 

70 

80 


INPUT - World ocean depth profile 
INPUT - Initial tide solution 
OUTPUT - Tide solution cotidal and corange 
values 

OUTPUT - Time series of tide heights over 
the span of integration for three 


INPUT ■ 
OUTPUT 


specified geographic location. 

^ 12 IS. 

' ax’ 3(j)’ ax’ a<i) 

- Tide solution 5 at the final integration 

time. 


Module II 30 

40 
50 


INPUT - Tide solution cotidal and corange 


values (from unit 50 of Module I) 
OUTPUT - Tide perturbation potential T' 


OUTPUT 


in cotidal and corange form, 
ap 3P 35. 

ax’ 3!j)’ ax’ 3(j) 




Table IX 


Vf^riable 

Narn e 

DOELT 

TTK 

BR 

SR 

AK 

K 


K2 

CT ; 
SIG 
G 
W 

m 

IFLGR 


MODULE I NAMELIST VARIABLES 


Function 


Integration time step in minutes. 

Integration start time in minutes when the restart 
option is used. 

Radius of the earth in meters. 

Friction coefficient, C . 


Default Value 

( 6 .) 

( 0 .) 

(6371 X 10^) 
(.003) 


Coefficient of lateral turbulent viscosity , (1 x 10 ) 

meters /sec. 


Dimensionless "general lunar coefficient" 


(.842 X 10“^) 


M 


moon 


2 M 


earth 


(^) 

y moon f 


x^here M denotes mass , the earth ’ s radius and 

R the mean lunar distance, 
moon 

Love number which represents the coefficient of 
additional potential due to deformation by a 
non-loading potential. 


( 0 .) 


(.1405 X 10“^) 


Amplitude factor for tidal constituent. (C of Eq (25)) (!•) 

Tidal constituent frequency in rad/sec. 

, , , _ 

Gravitational acceleration in meters/sec . 

Rotational rate of the earth in rad/sec. 

Grid size in degrees. 


Restart Flag 

0 - Start integration from t=0 from zero initial 

solution. 

1 - Read tidal solution at time TTN from I/O Unit 

40 and integrate over 1/4 tidal period; Form 
corange and co tidal solutions. 

2 - Read tidal solution at time TTN from I/O Unit 

40 and integrate over NITR tidal periods. 


(9.81) 

(.72722 E-4) 
(6^) 

( 0 ) 
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Table IX (Co at:.) 


MODULE I NAMELIST VARIABLES 


Variable 
Name 

IPERT 


^ Function Default Value 

Perturbation Flag (,U> 

0 - No perturbation potential F’ included in 

equations of motion. 

1 - Read |^, *|^ from I/O Unit 70, and 

include F’ potential in equations of motion. 


NOO Lower latitude limit in degrees for global inte- (-78®) 

gration. NOTE: The difference between -90° and 

NOO roust be an integer times the grid size ND. 

NEE Upper latitude limit in degrees for global inte- (84°) 

gration. NOTE: The difference between 90° and 

NEE must be an integer times the grid size ND. 


NITR 

The number of tidal constituent periods over vzhich 
the integration is to proceed. 

(1.) 

N1 

Grid point number (in compacted form) to be sampled 
for a time series analysis. 

(10 

N2 

Grid point number (in compacted form) to be sampled 
for a time series analysis. 

(2 0 

K3 

Grid point number (in compacted form) to be sampled 
for a time scries analysis. 

(3.) 

ITYPE 

Species of tidal component being integrated. 

2 - diurnal 

3 - semi-diurnal 

(30 

H2 

Love number which represents the coefficient of 
uplift response of the. solid earth to a non-loading 
potential. 

\ (O.) 


Table X 

MODULE II NAMELIST VARIABLES 


Variable 

Name 

Function 

Default Value 

BR 

Radius of the earth in meters (should be the 
same value used to generate the tidal solution 
from Module I. 

(6371x10^) 

NOO 

i' 

Lower latitude limit in degrees for global integration, (-78®) 
(Must be the same value used to generate the 
tidal solution from Module I. 

1 NEE 

Upper latitude limit in degrees for global 
integration. (Must be the same value used to 
generate the tidal solution from Module I. 

(84») 

! GFCNG 

Greens Function Flag 

(0.) 


0. - Use default Greens funtion 

(See Table III) 


1 ■ 

1. - Use input values of and U* functions 

to calculate Greens function. 


GFHP 

U' Input Flag (Valid ONLY IF GFING=1.) 

(1.) 


0, “ Set U’ to zero in the Greens Function 



1. - Input U’ as a function of angle (0®-180®) 

in increments of 1° [F0RMAT 8F10.3] in units 
U’ (q)*A*0*1O^2. 


GFKP 

Input Flag (Valid ONLY IF GFGNG = 1. ) 

(1.) 


0, - Set to zero in the Greens Function 

1, luput as a function of angle (0® - 180®) 

in increments of 1° [F0RMAT 8F10.3] in units 
$^(0) * A *0* 10l2. 


GFSG 

Self-Gravitation Term Flag 

(1.) 


0. - Set self-gravitation term to zero in the 

Greens Function. 



1. - Include self-gravitation term in the 

Greens Function. 
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5. 2 Prograni Output 


Ih. output fro. tho lutegrotiop algorid- ovut th. oco.a grid a.tworlc 

IS voluminous and difficult to analyze in raw for.. This naceasitatas plotted 
displays for visual analysis of th. solution. I. particular, the principal outputs 
are th- cotldal and corange solution contours and the tidal time series plots for 
selected geographic point, which are output to a plot tape for interface with an 
eatemal plotting device. In addition, the Integrator time t and the variables 
V, C and h+s are output to tape in their compacted ocean arrays at the end of 
elch integration period for program restart purposes. Additionally, the tidal 
time series is output on tape at the ternination of the integration. 

Printed plots are also produced for every plot output to the eac.rnal plotting 
device, in addition, the progra. prints th. total world function and the world 
grid of the ocean depth values for user inspection prior to calling the Integra 
non algorithm. At the end of each integration period the variables u. v and i 
are printed in their compacted ocean arrays, and then are printed in a full world 

grid array display. 
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